#include "ljforce.h"
#include "box.h"
#include "parameters.h"
#include <algorithm>

void LjForce::LjTerms(const double &r, double &attractive_term, double &repulsive_term)
{
    attractive_term = pow((sigma_ / r), 6);
    repulsive_term = pow(attractive_term, 2);
    return;
}
double LjForce::ForceFunc(const double &rsquare, const double &attractive_term, const double &repulsive_term)
{
    double fij = 4.0 * epsilon_ * (12.0 * repulsive_term - 6.0 * attractive_term) / rsquare;
    return fij;
}
double LjForce::PotentialFunc(const double &attractive_term, const double &repulsive_term)
{
    double potential_energy = 4.0 * epsilon_ * (repulsive_term - attractive_term);
    return potential_energy;
}
void LjForce::Distance(const std::vector<double> &ri, const std::vector<double> &rj, std::vector<double> &rij, double &r, double &r2)
{
    for (int k = 0; k < 3; k++)
    {
        rij[k] = ri[k] - rj[k];
    }
    if (gIsPeriodic) {
        PutInsideBox(rij, gBoxLength);
    }
    r2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
    r = sqrt(r2);
    return;
}

void LjForce::ComputeForce(MolSys &mol)
{
    int natoms = mol.natoms_;
    // initialize force
    for (int i = 0; i < natoms; i++)
    {
        mol.atoms_[i].SetForce({0, 0, 0});
    }
    std::vector<double> rij(3);
    double r = 0, r2 = 0;
    double attractive_term = 0;
    double repulsive_term = 0;
    double potential = 0, virial = 0;
    // cutoff radius, and its square
    double rcutoff = std::min(gBoxLength / 2.0, 2.5 * sigma_);
    double rcutoff2 = rcutoff * rcutoff;
    for (int i = 0; i < natoms - 1; i++)
    {
        for (int j = i + 1; j < natoms; j++)
        {
            Distance(mol.atoms_[i].coord_, mol.atoms_[j].coord_, rij, r, r2);
            if (r2 < rcutoff2)
            {
                LjTerms(r, attractive_term, repulsive_term);
                double fij = ForceFunc(r2, attractive_term, repulsive_term);
                for (int k = 0; k < 3; k++)
                {
                    mol.atoms_[i].force_[k] += fij * rij[k];
                    mol.atoms_[j].force_[k] -= fij * rij[k];
                }
                double energyij = PotentialFunc(attractive_term, repulsive_term);
                potential += energyij;
                virial += 4.0 * epsilon_ * (12.0 * repulsive_term - 6.0 * attractive_term);
            }
        }
    }
    mol.potential_energy_ = potential;
    return;
}